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GENERATION  OF  PSEUDO-RANDOM  NUMBERS 


I.  INTRODUCTION 


The  advent  of  high-speed  digital  computers  has  made  it  possible  to  use  simula¬ 
tion  techniques  incorporating  probabilistic  features.  These  simulations  are  generaUy 
referred  to  as  Monte  Carlo  simulations  and  are  resorted  to  whenever  the  systems  being 
studied  are  not  amenable  to  deterministic  analytical  methods  or  where  direct  experi¬ 
mentation  is  not  feasible.  An  integral  part  of  these  simulations  is  the  use  of  random 
numbers  having  a  certain  specified  distribution  characteristic  of  the  process  being 
studied.  Since  these  simulations  require  a  vast  amount  of  calculations,  speed  is  a 
vitally  important  factor.  It  was  soon  recognized  that  it  was  intractable  to  fill  the 
computer's  memory  with  large  aggregates  of  random  numbers.  Because  of  this, 
computer  algorithms  were  developed  which  allowed  random  numbers  to  be  generated 
on-line.  Since  the  generation  of  random  numbers  by  such  numerical  algorithms  is 
somewhat  a  contradiction  in  terms,  they  are  often  called  "pseudo -random"  numbers. 
There  are  essentially  two  problems  encountered  in  generating  such  random  numbers. 
One  is  that  the  generated  random  numbers  are  not  representative  of  the  desired 
distribution,  and  the  other  is  that  they  are  not  statistically  random,  i.e.,  that  there 
exist  correlations  in  the  generated  numbers.  The  latter  problem  is  the  more  serious 
one  as  evidenced  by  the  considerable  attention  it  has  been  given  in  the  literature. 

This  report  presents  several  random  number  generators  which  have  been  found 
particularly  useful  in  aerospace  engineering  applications.  APL  computer  programs  are 
also  listed  in  Appendix  B  for  most  of  the  generators. 


II.  UNIFORM  RANDOM  NUMBER  GENERATORS 


The  basic  element  of  all  Monte  Carlo  simulations  is  the  uniform  random  number 
generator.  Once  uniform  random  numbers  are  available,  all  other  desired  distributions 
can  be  obtained  either  by  use  of  the  probability  integral  transformation  or  by  applying 
some  known  relationship  between  the  desired  distribution  to  be  generated  and  the 
uniform  distribution,  1  One  of  the  early  methods  used  to  generate  uniform  random 
numbers  was  the  mid-square  method  originally  proposed  by  John  V.  Neumann,  In 
practice,  one  selects  an  arbitrary  K-digit  number,  squares  it,  and  then  selects  the  K 
middle  digits  as  the  new  random  number.  The  process  is  repeated  using  this  new 
random  number.  The  drawback  of  this  method  is  that  it  can  produce  a  zero  random 
number  at  unpredictable  times  and,  thus,  the  process  terminates.  Consequently, 
this  method  was  abandoned  quite  early  in  favor  of  the  so-called  congruential  method 
first  proposed  by  D,  H.  Lehmer  in  1949  [1].  Accordingly,  the  random  number 
generator  takes  the  form 

i;  Often  it  is  convenient  to  generate  a  random  number  from  a  specified  distribution 
by  employing  its  relationship  to  the  normal  distribution.  However,  one  must  first 
generate  the  normal  random  number (s)  as  a  function  of  uniform  random  numbers 
and  then  proceed  to  the  desired  distribution. 


(1) 


=  (a  X.  +  b)  Mod  M 
1+1  1 


where  the  multiplier  a,  increment  b,  and  modulus  M  are  integers.  The  starting  value 
X^  is  called  the  "seed"  of  the  random  number  generator.  Since  the  congruential 

relationship  (1)  is  cyclical,  the  sequence  of  random  numbers  will  repeat  after  a  cer¬ 
tain  period.  Generators  in  which  b  =  0  are  called  multiplicative;  otherwise,  they  are 
called  mixed.  The  statistical  behavior  of  the  generated  random  numbers  is  predomi¬ 
nantly  governed  by  the  choice  of  the  multiplier  a  and  the  modulus  M.  Therefore,  the 
most  widely  used  generators  are  of  the  multiplicative  type  (b  =  0) .  Many  empirical 
and  theoretical  tests  have  been  developed  to  assess  the  "goodness"  of  a  random 
number  generator.  One  of  the  most  popular  of  these  tests  today  is  the  lattice  test 
which  determines  the  lattice  structure  of  the  random  number  generator  by  comparing 
successive  n-tuples  (x. ,  •••>  ^^i+1’  ^i+2’  ^i+n^’  According 

to  this  test  an  acceptable  generator  can  be  obtained  by  selecting  the  constants  a,  b, 
and  M  to  achieve  a  nearly  hyper-cubic  lattice  structure,  i.e.  ,  the  ratio  of  cell  sides 
should  be  close  to  unity  [2].  In  practice,  n  is  usually  less  than  or  equal  to  five  [3]. 
Especially  useful  are  congruential  generators  for  which  the  modulus  M  is  a  prime 
number  and  hence,  are  called  prime  modulus  generators.  If  the  multiplier  a  is  selected 
to  be  a  primitive  root  modulo  M,  the  generated  random  number  sequence  attains  its 
maximum  period  P  =  M-1  and  all  possible  value  from  0  to  M-1  will  be  generated. 

Two  useful  uniform  random  number  generators  of  this  type  which  have  very 
satisfactory  lattice  structures  are: 


a  =  7^  =  16,807  ,  b  =  0  ,  M  =  2^^-l  ^  2,147,483,647 


(2) 


with  5-dimensional  cell  side  ratios  1 ;  7.  60 : 3.  39: 2. 09 : 1.  67 ,  and 


a  =  7^02,479  ivio(j(2^^-l)  ^  29,903,947  ,  b  =  0  ,  and 

M  =  2^^-l  =  2,147,483,647  ,  (3) 


with  5-dimensional  cell  side  ratios:  1:1.04:1.30:1.22:1.09  [4], 

III.  RANDOM  NUMBERS  FROM  CONTINUOUS  DISTRIBUTIONS 


Many  of  the  generators  for  continuous  distributions  are  obtained  by  a  direct 
application  of  the  probability  integral  transformation  [5].  For  a  given  uniform  random 
number  u  between  zero  and  one,  a  random  number  x  having  the  desired  distribution 
F(x)  is  obtained  by  solving  the  equation  u  =  F(x)  for  x.  Since  this  process  requires 

the  determination  of  the  inverse  cumulative  distribution  function  F  ^(x)  ,  its  prac¬ 
ticality  depends  upon  the  availability  of  explicit  expressions  or  convenient  approxima¬ 
tions  for  this  inverse  cumulative  distribution  function.  We  now  discuss  methods  on 
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how  to  generate  random  numbers  from  continuous  distributions  which  appear  frequently 
in  aerospace  engineering  simulations. 


A.  The  Normal  Distribution 

The  most  common  distribution  is  the  normal  or  Gaussian  distribution  with  den¬ 
sity  function  given  by 


f(x;ii  ,0^) 


/  2t\o^ 


-00  <  X,  y  <  “  and  a  >  0 


(4) 


Consequently,  a  variety  of  methods  have  been  devised  which  can  be  used  to  generate 
normal  random  numbers.  Of  these,  the  following  three  were  found  to  be  satisfactory 
and  practical  in  terms  of  accuracy  and  computer  run  time. 

1.  Hastings  Approximation 

This  method  invokes  the  probability  integral  transformation  in  a  slight  variation 
in  that  it  employs  the  complement  of  the  cumulative  distribution  Q(x)  =  l-F(x).  The 
reason  for  this  is  that  suitable  rational  approximations  for  Q(x)  were  derived  by 
Hastings  [6].  The  most  accurate  of  Hastings  approximations  is  given  as: 


X  =  t 


2 

C  +  C,  t  +  t 

O  1  2 

1  +  d^t  +  d2t^  +  dgt^ 


+  e(p) 


(5) 


where 


C  =  2.515517 
o 

=  0.802853 
C„  =  0.010328 
dj^  =  1.432788,  t  =  /- 21np 
d2  =  0.189269,  0  <  p  <  0.  5 
d„  =  0.001308,  and  ]  e(p)  |  <  4.  5  x  10 


X  is  the  desired  Normal  random  number  and  p  is  a  uniform  random  number. 

2.  Box-Mueller  or  "Polar”  Method 

This  method  generates  a  pair  of  normal  random  variables  using  a  pair  of  uni¬ 
form  random  numbers  as  follows: 
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Let  and  U2  be  independent  uniform  random  variables  and  define 
1/2 

=  (-21n  u^)  cos  2Tr  , 

(6) 

1  /2 

X2  =  (-21n  sin  2tj  U2 

Then  x^^  and  X2  are  two  independent  normal  random  variables  with  zero  mean  and  unit 
variance.  To  see  this,  we  establish  the  inverse  relationships 


(7) 


probability  density  function  of  x^,X2  is 
f  (x^,X2)  =  ^  =  f(x^)  f(x2)  ,  (8) 


^  2  ^  2, 
(x^  +  X2  ) 


u^  =  e 


and 


-1  ,  -1  /  ^2 
Uo  =  K-  tan  — 
2  2tt  I  x^ 


It  follows  then  that  the  joint 


and,  thus,  the  desired  conclusions,  including  the  independence  of  x.  and  x„,  are 
obtained  [7].  ^ 

3.  Central  Limit  Method 


Let  x^,  X2, 


x^  be  a  sequence  of  n  uniform  random  variables.  Then 


(9) 


will  be  distributed  asymptotically  as  a  normal  random  variable  with  zero  mean  and  unit 
variance.  For  n  =  12,  we  see  that  (9)  reduces  to 


y 


12 


6 


4 


The  exact  distribution  of  the  standardized  sum  of  n  independent  uniform  random 
numbers  can  be  easily  derived  using  moment  generating  functions.  Since  the  sum  x 
of  n  independent  uniform  variables  has  moment  generating  function 


(10) 


its  density  function  is  given  as 


f(x) 


1 

(n-1)! 


X]  (k)  uCx-k)]"^-^ 

k=0 


(11) 


where 


u(s)  = 


0 

1 


s  <  0 

s  >  0 


and  0  i  X  <  n 


For  the  standardized  random  variable  y^  in  equation  (9),  we  find  that 


(y) 


172 


n 


(n-1) 


(-1) 


k 


k=0 


+ 


n 

2 


where 


i/Sn 


i/5n 


(12) 


Density  functions  obtained  by  this  method  for  n  =  2,  4,  12,  20  are  compared 
with  the  normal  density  function  (dotted)  in  Figure  1.  Also,  a  comparison  of  the 
cumulative  distribution  function  for  n  =  12  with  the  cumulative  normal  distribution  is 
given  in  Table  1  to  four  decimal  places. 

The  agreement  between  the  two  distributions  is  very  good  except  in  the  tail 
areas.  But  a  comparison  of  random  numbers  generated  by  the  three  methods  revealed 
no  significant  statistical  differences  even  for  the  tail  areas. 

A  comparison  of  computer  CPU  run  times  for  each  of  these  three  methods  to 
generate  1000  normal  random  numbers  is  as  follows: 
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Method 

Time 

Hastings  Approximation 

Box-Mueller 

Central  Limit  Method 

1.13  sec 

1.02  sec 

1.08  sec 

Figure  1.  Normal  density  function  (dotted)  and  standardized 
sum  of  uniforms  (n  =  2,  4,  12,  20). 

TABLE  1.  CUMULATIVE  NORMAL  DISTRIBUTION 


X 

F(x) 

-4.0 

3.167  X  10'^ 

8.526  X  10'® 

-3.  5 

2.326  X  10'^ 

1.212  X  10'^ 

-3.0 

1.350  X  10'^ 

1.007  X  10"^ 

-2.  5 

6.210  X  10"^ 

5.579  X  10"^ 

-2.0 

2.275  X  10'^ 

2.227  X  10"^ 

-1.5 

6.681  X  10'^ 

6.  745  X  10'^ 

-1.0 

1.587  X  10'^ 

1.608  X  10"^ 

-0.  5 

3.085  X  10“^ 

3.106  X  10“^ 

0.0 

5.000  X  10'^ 

5.000  X  10"^ 

0.5 

6.915  X  10'^ 

6.894  X  10“^ 

1.0 

8.413  X  10"^ 

8.393  X  10"^ 

1.5 

9.332  X  10"^ 

9.326  X  10'^ 

2.0 

9.773  X  10'^ 

9.777  X  10"^ 

2.  5 

9.  938  X  10'^ 

9.944  X  10'^ 

3.0 

9.987  X  10’^ 

9.  990  X  10'^ 

3.  5 

9.998  X  10’^ 

9.999  X  10'^ 

4.0 

1,000 

1.000 

♦Comparison  of  the  cumulative  normal  distribution 
F(x)  with  the  cumulative  distribution  function  for 
the  central  limit  method  using-  equation  (13)  as  the 
density  function  with  n  =  12. 


B .  Log-Normal  Distribution 

It  is  often  claimed  that  the  log-normal  distribution  is  as  fundamental  as  the 
normal  distribution  and  may  be  thought  of  as  arising  from  the  combination  of  random 
terms  by  a  multiplicative  process.  The  log-normal  distribution  has  been  applied  in  a 
wide  variety  of  fields  including  social  sciences ,  physical  sciences ,  and  engineering  and 
its  density  function  is  given  by 


f(x;y,a)  =  (2tto2x2)'^/^  exp  ["-1/2  ( ""a"  )  ]  ’  0  <  x  ,  a  <  -  (13) 

I-  ^  -oo<y<oo^ 


Log-normal  random  numbers  x  may  be  generated  by  the  relationship  x  e  ,  where  y 
is  a  normal  random  number  obtained  by  methods  discussed  in  Section  III,  paragraph  A 


C.  Weibull  Distribution 

The  Weibull  distribution  is  perhaps  the  most  popular  distribution  at  the  present 
time  when  dealing  with  problems  of  reliability  and  material  fatigue.  Its  appeal  stems 
from  its  mathematical  tractability  and  for  this  reason  is  often  preferred  to  the  gamma 
distribution.  The  Weibull  density  function  is  given  by 


f(x;a,3)  =  a3  x^'^  exp(-a  x^)  ,  0  <  x  <  “  a,  g  >  0  .  (14) 

The  cumulative  distribution  function  F(x)  -  1  leads  immediately  to  the  inverse 

relationship 

^  -  \^lnu  ) 

as  the  desired  Weibull  random  generator. 


D.  Gamma  Distribution 

The  gamma  distribution  is  another  two  parameter  distri.bution  which  is  also  quite 
flexible  in  fitting  a  variety  of  random  processes.  It  finds  use  in  reliability  analysis, 
meteorology,  and  aerospace  engineering.  Its  density  function  is 


f(x)  = 


-ox  g- 1 

e  X 


(16) 


where  0  <  x  <  “>  and  a,g  >  0, 
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For  integer  values  of  3,  the  gamma  distribution  is  often  referred  to  as  the  Erlangian 
distribution  after  the  Danish  mathematician,  A.  K.  Erlang,  who  introduced  it  in  the 
theory  of  queues  and  Markov  processes  in  1917.  Random  numbers  following  the 
Erlangian  distribution  are  generated  by  the  formula 

3 

i=l 


E.  Exponential  Distribution 

The  exponential  distribution  appears  often  in  engineering  applications  because  of 
its  importance  in  reliability  theory  and  queueing  theory.  The  probability  density 
function  of  the  exponential  distribution  is 


f  (x ;  a)  =  a  e 


■ax 


(20) 


where  x,a  >  0.  In  this  case,  an  explicit  expression  for  the  inverse  cumulative  dis¬ 


tribution  function  is  easily  obtained  by  solving  the  equation  u  —  1  -  e  for  x, 
producing 


In  u 

a 


(21) 


as  the  desired  exponential  random  number  generator. 


F.  Chi-Square  Distribution 

2 

As  a  special  case  of  the  gamma  distribution,  the  x  -distribution  is  often  used 
as  a  measure  of  goodness  of  fit  of  a  specified  distribution  to  observed  frequencies. 

2 

For  discrete  variates,  x  provides  a  sensitive  test  of  departure  from  the  Poisson 
distribution.  The  Chi-Square  density  function  with  n  degrees  of  freedom  is  given  by 


X  n  1  x/ 2 

f(x;n)  =  - - -  ,  X  >  0  and  n  ^  1,2,3,...  .  (22) 

2  r  (n) 


Because  of  its  relationship  to  the  normal  distribution,  Chi-Square  random  numbers 
may  be  generated  by  taking  the  sum  of  n  squared  Normal  random  numbers. 
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G .  F  Distribution 


The  F  distribution  appears  extensively  in  statistical  hypothesis  testing  under 
normality  theory  and  has  the  density  function 


f(x;m  ,n) 


^n/2  -  1 

r(n/2)  r(m/2)  (m 


(23) 


where  x  >  0  and  m,n  =  1,2,...,  are  degrees  of  freedom.  In  practice,  it  is  generally 
found  to  be  most  convenient  to  make  use  of  the  F  distribution's  relationship  to  the 

2 

Chi-Square  distribution  in  order  to  generate  F  random  numbers.  That  is,  if  x  (™) 

and  x^(^)  two  independent  Chi-Square  random  variables  with  m  and  n  degrees  of 
freedom,  respectively,  then 


Y 


X^(m)ym 

X^(n)7n 


(24) 


follows  the  F- distribution  with  m,  n  degrees  of  freedom.  Thus, 


m 


i=m+l 


(25) 


is  the  desired  F  random  number  generator  with  m  and  n  degrees  of  freedom  and  the 
X.,  i  =  l,2,...,m+n  are  Normal  random  numbers. 


H.  Beta  Distribution 

To  generate  random  numbers  from  the  Beta  distribution  with  density  function 
given  by: 


f(x;a,  3) 


r  (g  +  3)  a-1 

r(a)  1(3) 


(1  -  X) 


3-1 


0  <  X  <  1 


0  <  a,  3  <  “ ,  (26) 


it  is  most  convenient  to  make  use  of  a  relationship  between  the  Beta  distribution  and 
the  F- distribution.  That  is,  if  X  is  a  random  variable  following  the  F  distribution 

with  m  and  n  degrees  of  freedom,  then  Y  =  ^  ^  follows  the  Beta  distribution. 

Consequently, 
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Y 


2m 


X. 

1 


2 


2n+2m 


i=2m+l 


(27) 


is  the  desired  Beta  random  numbers  generator,  where  the  Xds  are  standard  Normal 
random  numbers  realized  by  methods  in  Section  III,  paragraph  A. 

IV.  RANDOM  NUMBERS  FROM  DISCRETE  DISTRIBUTIONS 


Generation  of  random  numbers  from  a  discrete  distribution  is  handled  in  a 
manner  analogous  to  that  of  a  continuous  distribution.  Given  a  uniform  random  number 
u  and  cumulative  discrete  distribution  F(x) ,  the  least  value  of  x  for  which  F(x)  >u  is 
sought.  This  X  is  the  desired  random  number  having  discrete  density  function  f(x). 


A.  Binomial  Distribution 
The  cumulative  binominal  distribution  is  given  by 


X 

F(x;n,p)  = 

k=0 

For  a  given  uniform  number  u,  one  may  successively  evaluate  the  upper  limit  x  until 
the  minimum  value  of  x  is  found  for  which 


(S) 


k 


(1  -  P) 


x-k 


(28) 


X 

i:( 

k=0 


k  .n-k  ^ 

p  (1  -  p)  A  u 


(29) 


An  alternate  method  is  based  on  the  sum  of  n  independent  Bernoulli  random  variables. 

By  this  method,  Bernoulli  trials,  each  having  probability  of  success  p,  are  simulated 
with  each  Bernoulli  trial  being  assigned  a  one  or  zero  (depending  on  success  or  failure). 


2.  Because  factorials,  exponentials,  and  powers  frequently  occur  in  probability 
density  functions,  care  must  be  taken  to  avoid  computer  overflows /underflows 
when  computing  individual  components  of  the  density  function.  Recursive  rela¬ 
tionships  are  quite  useful  in  dealing  with  these  types  of  problems,  but  computer 
programs  employing  recursive  relationships  are,  in  general,  much  slower. 
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Then  the  siim  of  these  n  Bernoulli  random  numbers  will  be  a  random  number  from  a 
Binomial  population  with  parameters  n  and  p.  However,  it  was  found  that  this  process 
is  relatively  slow. 


B .  Poisson  Distribution 

Poisson  random  numbers  x  may  be  generated  from  cumulative  distribution 
function 


F(x;A) 


^  ,k  -A 

\  ^  A  e 

Lj  k! 

k=0 


A  >  0 


(30) 


by  using  the  same  technique  which  was  used  for  the  binominal  distribution.  An 
alternate  method  for  obtaining  Poisson  random  numbers  is  by  generating  uniform 
random  numbers  u^  until  the  inequality 


k+1 

TT 

i=l 


u.  <  e 
1 


(31) 


is  satisfied,  which  gives  k  as  the  desired  Poisson  random  number.  However,  this 
method  was  found  to  be  relatively  slow. 


C.  Neyman  Type-A  and  Thomas  Distributions 

These  two  distributions  belong  to  the  category  of  cluster  (self- exciting)  point 
processes  and  have  found  application  in  aerospace  engineering,  ecology,  reliability, 
and  forestry.  Cluster  processes  are  characterized  by  a  primary  (mother)  process 
which  generates  at  each  point  secondary  (daughter)  events.  When  only  daughter 
events  appear  in  the  final  process  and  when  primary  and  secondary  distributions  are 
both  Poisson,  the  resulting  distribution  is  known  as  the  Neyman  type- A  counting 
distribution  with  probability  function  [  8] 


p(n;a,3)  =  ^ 


.  .n  -am 
(am)  e 

n! 


m! 


m=0 


n  =  0,1,2,. . . 
and  a, 3  >  0 


(32) 


The  parameter  3  represents  the  rate  at  which  primary  events  occur  and  a  is  the 
average  number  of  secondary  counts  per  primary.  The  mean  and  variance  is  a3  and 
(1  +  a)a3,  respectively.  The  Thomas  distribution  is  similar,  except  that  primary 
events  are  counted  in  the  final  process.  The  probability  function  is  given  by 
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n-m 


n  =  0,  1 ,  2,  . . . 


(33) 


n 


p(n;a,  3) 


E 

m=0 


(am) 


„m  -  3 
-  am  3  e 


(n-m)! 


m; 


2 

with  mean  (1  +  a)  3  and  variance  (1  +  3a  +  a  )3. 

APL  programs  are  provided  in  Appendix  B  which  generate  random  numbers  from 
these  distributions. 


V.  CONCLUSION 


Because  large  scale  simulations  often  require  a  vast  number  of  random  numbers 
from  various  distributions,  emphasis  should  be  placed  on  speed  and  accuracy  of  the 
methods  used.  We  have  found  the  congruential  uniform  generators  presented  in 
this  report  to  be  both  suitable  and  practical.  However,  not  all  possibilities  have  been 
exhausted  and  the  literature  is  found  to  be  extensive  in  this  area.  Random  numbers 
can  be  generated  from  nonuniform  distributions  by  finding  the  inverse  to  the  cumula¬ 
tive  distribution  functions  in  accordance  with  the  probability  integral  transformation. 
However,  when  no  convenient  inverse  exists,  numerical  approximation  methods  or 
statistical  relationships  may  be  used  to  generate  the  desired  random  numbers.  Again, 
speed  and  accuracy  should  be  key  factors  in  choosing  between  candidate  methods. 

Some  useful  approximations  to  the  cumulative  normal  distribution  are  given  in  Appendix 
A  and  APL  programs  for  most  of  the  generators  are  listed  in  Appendix  B. 
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APPENDIX  A 


APPROXIMATION  TO  THE  CUMULATIVE  NORMAL  DISTRIBUTION 


Because  of  the  extensive  application  of  the  Normal  distribution  in  aerospace 
simulation  studies,  we  include  here  three  methods  for  approximating  the  cumulative 
standard  Normal  distribution.  Other  useful  formulae  may  be  found  in  Reference  9. 
Denote  the  standard  normal  density  function  by 


f(x)  =  (2tt) 


-1/2 


12 


and  its  cumulative  distribution  function  by  F(x).  Then  F(x)  may  be  represented  by 
each  of  the  following: 


1 

2 


oo 


E 

n=0 


(-1) 


n  2n+l 

X 


2%!  (2n+l) 


-OO  <  X  <  OO 


(A-1) 


If  If  ^  " 

,  r  2  2  „  2  „  2  .  2 

1  .  X  X  2x  3x  4x 

2  +t(,x)  ••• 


X  >  0 


(A- 2) 

(A- 3) 


The  last  two  formulae  are  called  continued  fractions,  where  the  conventional  notation 


bi+  b2+ 


bi+a2 


is  used  to  conserve  space.  Since  each  formula  is  exact,  we  investigate  their  rates  of 
convergence  to  F(x)  for  computational  convenience. 

Equation  A-2  converges  most  rapidly  for  x  >  3  while  A-1  and  A- 3  are  preferred 
for  0  <  X  <  3. 
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APPENDIX  B 


APL  PROGRAMS 


'V  R^N  CDFI  CDF 

[  1  ]  R<r\0 

[2]  L:R^R,+/CDF<URF  1 

[3]  -^{0<N^N-1) /L 

[  4  ]  -»-0 

[5]  ft  RETURNS  N  RANDOM  NUMBERS  FROM  TEE 

[6]  fi  DISCRETE  DISTRIBUTION  FUNCTION 

[7]  ft  SPECIFIED  BY  CDF. 

V 

V  R^CTRALIM  N 

[1]  R*-  6+(  +  /  l  +  ?Npl001)ilOOO;N-^N  ,12 

[2]  ^0 

[3]  ft  THIS  FUNCTION  GENERATES  NORMAL 

[4]  ft  RANDOM  NUMBERS  ACCORDING  TO 

[5]  ft  TEE  CENTRAL  LIMIT  METHOD. 

V 

V  P-f-X  CUMNORM  X-,C\S',Y  \A',MiT  iZ%S0 

[1]  P«^(O2)*"0  . 

[2]  L0:T^X^X^Z>(  l*lZ-*-l+K-M-*-M+l 

[3]  A^(2^Z)+  1+TiA 

[4]  5-<-5  +  S0**-S0x-JxA:x  (  l  +  2xM)  ^^^x2xl  +  2xil^ 

[5]  Y<-X+ZrY 

[6]  ■^iM<K)/L0 

[7]  P^O  0  0 

[83  P[l]^l-Px(*-IxJff2)^J 

[9]  P[2]^0.5+(7x(*-Jfx:5r^2)xZ^>l 

[10]  P[3]^0.5  +  (7x5 

[11]  -»-0 

[12]  0  THREE  FORMULAE  FROM  THE  HANDBOOK  OF  MATHEMATICAL 

[13]  ft  FUNCTIONS ,  EDITED  BY  ABROMOWITZ ,  1970  ,  ARE  PROVIDED. 

[14]  ft  THE  FIRST  TWO  ARE  CONTINUED  FRACTIONS  AND  THE  THIRD  IS 

[15]  ft  A  POWER  SERIES.  P[l]  (EQ.  26.2.14,  PAGE  932  ,  GIVES 

[16]  0  P(X),  X>0.  P[2]  ALSO  GIVES  P(J)  AND  IS  EQ .  26.2.15. 

[17]  ft  P[3]  IS  EQ.  26.2.10.  [l]  CONVERGES  MOST  RAPIDLY  FOR 

[18]  ft  X>3  WHILE  P[2]  AND  P[3]  ARE  PREFERRED  FOR  0<X<3. 

[19]  ft  PROGRAMMED  BY  L.  HOWELL,  19  MAY  82. 

V 

V  E^PARAM  GENLOGNORM  N ',MU ',V AR 

[1]  VAR^PARAMi2^',MU^PARAMll'\ 

[2]  E^*MU+{VAR*0 .3)^N0RM  N 

[  3  ]  ->-0 

[4]  ft  GENERATES  N  LOG-NORMAL  RANDOM 

[5]  0  NUMBERS  WITH  PARAMETERS  MU  AND 

[6]  ft  VAR  (MU  AND  VAR  ARE  THE  MEAN 

[7]  ft  AND  VARIANCE  OF  THE  DEFINING  NORMAL 

[8]  fl  DISTRIBUTION) . 


V 


V  R^HASTING  N iC iDiUOiSU iDE ;T 

[1]  C’'^2. 515517  0. 802853  0. 010328  ;Z?-^1  1.432788  0 . 1 89269  0 . 001 308 

[2]  l/0^-"0.5001+(.’^?pl000)*1000 

[3]  T'^(-2x®  |i/0)*0.5 

[4]  NV->-Gbil+T^Cl2:\+TKCldJ 

[5]  Z>E'^Z)C1]+2’xZ?C2]+2’xZ?[3]+2’xZ?[4] 

[6]  R^(^U0)^(T-NU^DE) 

[7]  -t-O 

[8]  ft  GENERATES  NORMAL  RANDOM  NUMBERS 

[9]  ft  USING  BASTING'S  RATIONAL  APPROXIMATIONS 

V 


V  CDS-i-I/  NELYMM  A^ZidjMlX lY^N  ;NPiRiS 
Cl]  NP^*-V^l-*-A 

[2]  N^O 

[3]  Ll:N^N+l;S<-\0iM->-0 

[4]  L2:X'>-l;M->-l+M 

[5]  II-^l 

[6]  L3:X<-X^AyMTN+l-II 

[7]  ->-(N^II^II+l) /L3 

[8]  Ar-*-Jx*-^xAf 

[9] 

[10]  L^iiY-^-I^WiM+l-J 

[11]  HM^J-^J+D/Lii^ 

[12]  Y-^Y^*-W 

[13]  S^S.X^Y 

[14]  -►(  {R<0)y(R>lE'l2)viM^5)) /L2iR-*-/~2^S 

[15]  NP^NP,+/S 

[16]  -►(0.999>  +  /il?P)/Ll 

[17]  CDF-*-+\NP 

[18]  -»-0 

[19]  ft  COMPUTES  THE  CUMULATIVE  NEYMAN-TYPE-A 

[20]  ft  PROBABILITY  DISTRIBUTION  FUNCTION. 

[21]  ft  15  AVERAGE  NO.  OF  PRIMARY  (MOTHER)  EVENTS 

[22]  ft  15  AVERAGE  NO.  OF  SECONDARY  (DAUGHTER)  EVENTS  PER  PRIMARY 

[23]  ft  PROGRAMMED  BY  LEONARD  HOWELL^  AUGUST  4.  1982. 

V 

V  R^NORM  5 

[  1]  i?-<-"6  +  0.0001x  +  /?5pl0000;5-«-5.-12. 

V 
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7  PDF^PDFFRCDF  CDFiY 
[1]  PDF-<-CDF-Y  ;Y-<-0  1<^CDF 

V 

V  CDF^POISSON  STiPOiPiK 

[1]  CDF^PO->-*-ST 

[2]  K^P^l 

[3]  LOzP^P^STiK 

[4]  CDF-^CDF  ,i~‘i-^CDF)+P^PO 

[5]  K^K+1 

[6]  ->-iO.S999>~liCDF) /LO 

[7]  -►O 

[8]  ft  RETURNS  THE  CUMULATIVE  POISSON 

[9]  ft  DISTRIBUTION  FUNCTION  WITH 

[10]  ft  PARAMETER  ST. 

V 

V  R^MU  POISl  KiF;X 

[1]  R^(*-MU)^{MU*X)i\X‘>-0,\K 

[2]  F^+\R 

[3]  ->-0 

[4]  fi  CUM.  VIST.  FUNCT.  FOR  POISSON 

V 

V  R^MU  P0IS2  KiFl\P\X 

[1]  X‘<-OiP-^R->-(*-MU)  i 

[2]  L:P*-MU^P^X^X+1 

[3]  R^R,P 

[4]  -*-iK*X)/L 

[5]  F1^+\R 

V 


V  R^POLAR  N 

[1]  i?-H(?i?)^i?-(2.rx/A?.0.5)pl07  3741824 

[2]  R-*-Nq{1  2o.O(02)x”1  0+i?)x(  2  20.X01  0  +  J?)*0.5 

[3]  -*-0 

[4]  ft  GENERATES  RANDOM  NUMBERS  USING 

[5]  ft  TEE  BOX-MU  ELLER  METHOD,  OFTEN 

[6]  ft  REFERRED  TO  AS  THE  POLAR  METHOD. 

V 

V  R^NP  RBINOl  NiA;IiUiK 

[1] 

[2]  LO  :  i4-^0 

[3]  L:f/-<-Cl  +  ?1001)^1000 

[4]  K^U^NPL2l 

[  5  ]  A'*~A+K 

[6]  ->-{NPL11^I*I+1)/L 

[  7 ]  R^R . A 

[8]  -*-iO*N->-N-l)/LO 

[9]  -►O 

V 
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V  RBIN02  N iFiMiP;PFiU iXiY 

[1]  X^0,iM;M*-NPLl']iP^NPL2l 

[2]  PF^(.X\M)^(P*X)xil-P)*M-X 

[3]  F^+\PF 

[4]  i/'S-Cl  +  ?^?pl001)Tl000 

[5]  R^+/U°.>F 

[6]  -i-O 

[7]  fl  GENERATES  BINOMIAL  RANDOM  NUMBERS 

[8]  n  USING  THE  INVERSE  METHOD. 

[9]  R  ffP[l]  IS  THE  SAMPLE  SIZE  AND 

[10]  fi  NPL2)  IS  THE  PROBABILITY 

V 

V  R^NP  RBIN03  NiM\P\U 

[1]  ;P-<-fl?P[2] 

[2]  l+?(ff,M)pl001)vl000 

[3]  R^+/U>P 

[4]  ->0 

[5]  fl  GENERATES  BINOMIAL  RANDOM  NUMBERS 

[6]  fi  >}S  THE  SUM  OF  BERNOULLI  TRIALS. 

[7]  R  fP[l]  IS  THE  SAMPLE  SIZE  AND 

[8]  R  ffP[2]  IS  THE  PROBABILITY . 

V 


7  R^MU  RP0IS3  NiA;IiK;U 

[1]  R<-\0;I^1 

[2]  LO:A*-liK^Q 

[3]  P:P^(~l+?1001)f 1000 

[4]  A^AxU 

[5]  K^K+1 

[6]  ->-(A>(*-MU)) /L 

[7]  R-^R,Cl+K) 

[8]  ^{N>I^I+1) /LO 

[9]  -»-0 

[10]  fl  GENERATES  POISSON  RANDOM  NUMBERS 

[11]  R  SI  TAKING  THE  PRODUCT  OF  UNFORM 

[12]  R  RANDOM  NUMBERS  AND  CHECKING  FOR 

[13]  R  FOR  THE  INEQUALITY  EXPi-MEAN) 

V 

V  R^MU  RPOIS^  N',X\AiPiF-,U‘,M 

[1]  X^O\A^P<-{*-MU) 

[2]  L-.A^MU^AiX-^X^l 
[  3  ]  P^P ,  A 

[4]  ■^<.X*\M*-MU+^y-MU*0 .3) /L 

[  5  ]  P^P ,  A 

[6]  F^+\P 

[7]  l  +  ?i7pl001)^1000 

[8]  R^+/Uo.>F 

[9]  >0 

[10]  R  USES  THE  INVERSE  CUMULATIVE 

[11]  R  METHOD  FOR  DISCRETE  RANDOM 

[12]  R  NUMBERS 


^STDUNFSUMlWll 

V  STDUNFSUM  N iK ;X i S ; Z ;PRT lY lUX ; SIG iXl lU 
Cl]  UX^NT2iSIG-<-(Nil2)*0 .5 

[2]  X^-UXiSIGiK^~l  +  iiN+l)  iY^\OiPRT<rxO 

[3]  L:  Z^Ul-K)xU^Xl>KiXl^UX+SIG^X 

[4]  S^SIG>^(-/ {K[N)^Z*N-l)iliN-l) 

[5]  Y^Y,S 

[6]  PRT^PRT,X 

[7]  ;!r-*-J+0.0  5 

[8]  ->-{X<UXiSIG) /L 

[9]  0  DRAW  Y  FS  PRT 
CIO]  ^-0 

Cll]  COMPUTES  THE  PROBABILITY  DENSITY 
Cl2]  Pi  FUNCTION  OF  TEE  STANDARDIZED 
Cl3]  fl  SUM  OF  N  INDEPENDENT  UNIFORM 
Cl4]  PI  RANDOM  VARIABLES  BY  TAKING  THE 
Cl5]  PI  INVERSE  OF  ITS  MOMENT  GEN- 
Cl6]  ft  ERATING  FUNCTION.  PROGRAMMED 
Cl7]  PI  BY  LEONARD  HOWELL  08/04/82. 

V 


VSTUDENTiaiV 

V  T->-V  STUDENT  P ;  C ;  Z?  ;i? ; ZP ;  G1 ; (72  ;  G3  ;  G4  ;  G ;  7P ;  01 ;  P2 

Cl]  C-«-2. 515517  0.802853  0 . 0 1  0  3  2  8  ;  0-«-l  1.432788  0.189269  0.001308 

C2]  Pf-C2x«P)*0  .  5  ;7P-^7*"l  ~2  "3  ”4 

C3]  XP^R-D\iD2iDl-r-^/Cy~R*0  1  2\D2^+/D^R*0  12  3 

C4]  G1«^(+/ZP*3  l)f4 

C5]  G2-^(  +  /(ZP*5  3  px5  16  3)^96 

C6]  G3^(+/3  19  17  15xZP*7  5  3  1)^384 

C7]  G4^(+/79  776  1482  “l920  ■945xZP*9  753  1)^92160 

C8]  G^G1,G2,G3,G4 

C9]  2’'<-ZP+  +  /GxI7P 

CIO]  -»-0 

Cll]  PI  THIS  FUNCTION  APPROXIMATES  THE  INVERSE  STUDENT -T 
Cl2]  PI  DISTRIBUTION  AND  IS  USEFUL  IN  AUTOMATIC  LOOK-UP  OF 
Cl3]  PI  CRITICAL  VALUES  IN  HYPOTHESIS  TESTING.  FORMULA 
Cl4]  PI  26.7.5,  P^GP  949  OF  ABRAMOWITZ'  S 

C15]  PI  MZS.S.MiTlQ.AL_FUNQ.ZlQM  IS  USED  IN  CONJUNCTION  WITH 
Cl6]  ft  FORMuTa  26.2.23  {APPROXIMATES  THE  INVERSE  NORMAL 
Cl7]  PI  DISTRIBUTION)  y  PAGE  933.  PROGRAMMED  BY  L.  HOWELL  y 
Cl8]  PI  JUNE  23  ,  19  82  . 

V 
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lTEOMASi{]'\^ 

V  CDF^TEOMAS  I NP  iM -,0  ;  E ;  ED  ;  I  ’,J  iP  iT  IJ  iJl 

[1]  M^INPlllPi  AV  NO  OF  MOTEER  EVENT Si.PER  MS) 

[2]  D^INPL2lPi  AV  NO  OF  DAUGETER  EVENTS  (COUNTS) 

[3] 

[4]  T^ELll 

[5]  CDF^Elll  ,Ell^+EL2^xED-^*-D 

[6]  1^2 

[7]  LliP^J^l 

[8]  E^E,(  H-£')xWfJ 

[9]  L2:P^P^DiJ 

[10]  -^(I>J^J  +  1) /L2 

[11]  P^P-x-ED 

[12]  2’^Pxfi’[2] 

[13] 

[14]  L3  :P-^Px  (  (JliJ)*IJ )  xIJxED^DxJ  1-<-J  +  1 

[15]  T^T+PxElJ+2l 

[16]  ^(I>JW  +  ll/L3 

[17]  CDF^CDFA  liCDF)+T 

[18]  I^J+1 

[19]  ^-(  0.9999>~1+C£>F)/L1 

[20]  ->0 

[21]  PI  RETURNS  TEE  CUMULATIVE  TEOMAS  DISTRIBUTION 

V 


SIUNFLDIV 

V  U^UNF  S-,N 

[1]  y-«-i  0 

[  2 ]  L: U^U , UNFDxUNFX^UNFM \ UNFRxUNFX 

[3]  -^(0<N^N-1)./L 

[4]  U*-SpU 

[5]  ->0 

[6]  fl  STANDARD  UNIFORM  NUMBER  GENERATOR 

V 


Vi/fffIC[D]V 
V  UNFIC  SEED 

[1]  14748  3647 

[2]  i/il?Fi?'<-299  03947 

[3]  UNFX^SEED 

[4]  UNFD^iUNFM 

[5]  ^0 

[6]  PI  INITIALIZES  PARAMETERS  FOR  TEE 

[7]  0  STANDARD  UNIFORM  GENERATOR. 


V 


S/UNFSUMim^ 

V  UNFSUM_N;K;X;SiZ;UiPRTiY 

[1]  X^O;K^  l+\iN+l) ;Y^\0;PRT^i0 

[2]  L:Z^{X-K)^U^X>K 

[3]  f 

[4]  Y^Y,S 

[5]  PRT^PRT ,X 

[6]  X^X+0.1 

[7]  -y{X<N)/L 

[8]  0  DRAW  Y  VS  PRT-Ril 

[9]  R  COMPUTES  INVERSE  OF  THE  MOMENT 
CIO]  R  GENERATING  FUNCTION  OF  THE  SUM  OF 
[11]  R  INDEPENDENT  UNIFORM  VARIABLES 

V 


VUNIFORMWV 
V  R^M  UNIFORM  N\I\UiAiS 
Cl]  i  0  ;  I-<-l ;  U-^SEED  ;  >3^1 6  8  0  7  ;  MOD^~  1  +  2*31 

C2]  L:U^MOD\A^U 
C3]  S^IMxUtMOD 
C  4 ]  R^R , S 
C5]  -^(N>I^I+1)  /L 

[6]  +-0 

C7]  R  STANDARD  UNIFORM  NUMBER  GENERATOR. 


V 
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